/*!  \file  BiMRRR.h
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      09/10/2013 02:23:49 PM
 *   \copyright GNU Public License.
 */
#pragma once

#pragma once
#include "BiSVD.h"
//BISVD(BiMRRR);// Bidiagonal svd with Divide and conquer algorithm
class BiMRRR;
class GK;
class Interval;
class GKBisector;
#define BIMRRR_EXTRACT_PARAMS const double * RS z, double * RS u,\
 double * RS v, const double norm2, const int n, const int * sig
typedef void (BiMRRR::*extractFunc)(BIMRRR_EXTRACT_PARAMS) const;

class BiMRRR : public BiSVD { 
public: 
	BiMRRR();
	BiMRRR(const char * name);
	virtual ~BiMRRR(){};
	virtual void compute(const Matrix & B, Matrix & U2, Matrix & Sigma, Matrix & V2T);
	int getN() const { return n; };
	void setForComputing(const Matrix & B, Matrix & U, Matrix & S, Matrix & V);
	void doneComputing();

	// gap tolerance parameter
	double getGapTol      () const { return gapTol      ;}	
	// relative accuracy for B^TB bisector
	double getBisRelAcc   () const { return bisRelAcc   ;} 
	// relative accuracy, before computing singletons
	double getRefineRelAcc() const { return refineRelAcc;} 
	GK         * getNat   () const { return nat         ;}
	GKBisector * getBis   () const { return bis         ;}
	Interval   * getAll   () const { return all         ;}

	void computeSingleton(const GK & parrent, int index) const ;
	void computeCluster  (const GK & parrent, const int clusterStart, const int clusterStop) const;
	void computeGoodChild(const GK & parent, const int clusterStart, 
					const int clusterStop, GK & child, double & shift) const;
	void checkBorders(const GK & gk, const int i) const;
	void recallculateBorders(const GK & child, const int clusterStart, 
					const int clusterStop, const double shift) const;
	void computeEnvelopeVector(const GK & left, const GK & right, double * RS s2) const;

private:
	void setDefaultSettings();
	void extractSingularVectors(BIMRRR_EXTRACT_PARAMS) const {
		(this->*extractF)(z, u, v, norm2, n, sig);
	}
	void extractNoGap(BIMRRR_EXTRACT_PARAMS) const;
	void extractVGap (BIMRRR_EXTRACT_PARAMS) const;
	void extractUGap (BIMRRR_EXTRACT_PARAMS) const;
	void extractUVGap(BIMRRR_EXTRACT_PARAMS) const;
	extractFunc extractF;
	void setExtractFunctionPointer() {
		extractF = incU == 1 ?
			(incV == 1 ? &BiMRRR::extractNoGap : &BiMRRR::extractVGap) :
			(incV == 1 ? &BiMRRR::extractUGap  : &BiMRRR::extractUVGap);
	}

	bool checkGapTol(const double l1, const double l2) const {
		return fabs(l1 - l2) / fmax(l1, l2) >= gapTol;
	}


	Interval * all;
	double   * RS u;
	double   * RS v;
	double   * RS s;
	GK         * nat;
	GKBisector * bis;
	double       gapTol;	// gap tolerance parameter
	double    bisRelAcc; // relative accuracy for B^TB bisector
	double refineRelAcc; // relative accuracy, before computing singletons
	int n              ; // size of the problem
	int incU, incV, jumpU, jumpV;
	int maxIter;
};



class BiMRRRPar : public BiMRRR { 
public: 
	BiMRRRPar(): BiMRRR("BiMRRR parallel") {canBeParallel = true;};
	virtual bool setNumberOfThreads(const int nThreads_){ 
	//	nThreads = canBeParallel ? nThreads_ : 1;
		nThreads = 1;
		return canBeParallel;
	};
	virtual ~BiMRRRPar(){};
	virtual void compute(const Matrix & B, Matrix & U, Matrix & Sigma, Matrix & V);
};


